Effects of small surface tension in Hele-Shaw multifinger dynamics: 
an analytical and numerical study 
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We study the singular effects of vanishingly small surface tension on the dynamics of finger com- 
petition in the Saff man- Taylor problem, using the asymptotic techniques described in [S. Tanveer, 
Phil. Trans. R. Soc. Lond. A 343, 155 (1993)] and [M. Siegel, and S. Tanveer, Phys. Rev. Lett. 
76, 419 (1996)] as well as direct numerical computation, following the numerical scheme of [T. Hou, 
J. Lowengrub, and M. Shelley, J. Comp. Phys. 114 , 312 (1994)]. We demonstrate the dramatic 
effects of small surface tension on the late time evolution of two-finger configurations with respect 
to exact (non-singular) zero-surface tension solutions. The effect is present even when the relevant 
zero surface tension solution has asymptotic behavior consistent with selection theory. Such singular 
effects therefore cannot be traced back to steady state selection theory, and imply a drastic global 
change in the structure of phase-space flow. They can be interpreted in the framework of a recently 
introduced dynamical solvability scenario according to which surface tension unfolds the structually 
unstable flow, restoring the hyperbolicity of multifinger fixed points. 

PACS numbers: 47.54.+r, 47.20.Ma, 47.20.Ky, 47.20.Hw 
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I. INTRODUCTION 

The displacement of a viscous fluid by a less-viscous 
one in a Hele-Shaw cell, the so-called Saffman- Taylor 
problem jjj, |2|, |], [|, ||, is a prototypical pattern forma- 
tion problem. Since the seminal work of Saffman and 
Taylor [Q a considerable effort has been aimed at under- 
standing both steady and unsteady interfacial patterns 
formed during this flow. The Saffman- Taylor problem is 
the simplest member of a wide class of interfacial pat- 
tern formation problems such as free dendritic growth, 
directional solidification, or chemical electro-deposition 
H 0, ||. As such, a theoretical understanding of Hele- 
Shaw flow may help elucidate generic behavior common 
to many pattern forming systems. Despite its relatively 
simple formulation and the large amount of work devoted 
to it, however, several aspects of interfacial dynamics in 
Hele-Shaw flow are still poorly understood, in particular 
concerning the highly nonlinear and nonlocal dynamics 
of finger competition. 

One of the main reasons for the wide interest in Hele- 
Shaw flow, at least from a mathematical point of view, 
is that explicit time-dependent solutions can be found in 
the case of zero surface tension [§, 0, |§ . However, 
it is also known that the zero surface tension Saffman- 
Taylor (ST) problem is ill-posed as an initial value prob- 
lem |13| and finite-time singularities appear frequently 
jl4j]. Nevertheless, rather large classes of zero surface 
tension solutions have been found which exhibit the va- 
riety of morphologies observed both in experiments and 
numerical simulations. Then, the question that naturally 
arises is to what extent smooth (nonsingular) zero surface 
tension solutions reproduce the dynamics of the physical 
problem with finite surface tension, in particular in the 
limit of vanishing dimensionless surface tension, B — > 0. 

It is well known that surface tension is a singular per- 



turbation to the zero surface tension problem This 
singular character manifests dramatically in the classical 
selection problem posed by Saffman and Taylor jlj] and 
only solved three decades later (l^, [if], [l7], [l8| , where an 
arbitrarily small surface tension selects out a single, sta- 
ble solution from the continuum of steady single-finger 
B = solutions. Another manifestation of the singular- 
nature of surface tension which is directly relevant to the 
present work is its effect on the dynamics. Siegel, Tanveer 
and Dai |2(J showed that interfacial evolution for the 
regularized problem (i.e., vanishingly small B) may differ 
significantly from that for the B = problem in order 
one time. Then, smooth time dependent solutions of the 
B = case do not coincide, in general, with the limiting 
solutions for B — * 0. Accordingly, the study of finite B 
dynamics encounters considerable difficulties since B = 
solutions cannot be naively used as a starting point for 
the study of the problem with finite B. 

The physical content of exact zero surface tension so- 
lutions with pole-like singularities has been recently ad- 
dressed in Refs. || |T], |2^] using a dynamical systems 
approach. Through a detailed study it has been shown 
that the exact zero-surface tension phase flow, consid- 
ered in a global sense, is structurally unstable. In other 
words, the zero surface tension phase dynamics are not 
topologically equivalent to the phase space flow of the 
physical problem, regularized by surface tension. Indeed, 
the zero surface tension phase flow omits the necessary 
saddle-point structure of multifinger fixed points, which 
is crucial to the physical finger competition process [^2| . 
A natural extension of the well known solvability mech- 
anism (first applied to 'select' a finger of width 1/2 out 
of a continuum of solutions in the single finger case) was 
proposed for multi-finger solutions in |22j ; this helps clar- 
ify how the introduction of surface tension modifies the 
global phase space structure of the flow and restores the 
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hyperbolicity of multi-finger fixed points. 

however, was qualitative 



The approach of Ref. 22 



in nature and could not quantify the extent to which 
zero surface tension trajectories might resemble the evo- 
lution with small surface tension. In particular it was 
recognized that, while some trajectories appear to be 
qualitatively correct for infinite time, others may have 
a dramatically different evolution. In particular, adding 
an infinitcssimal surface tension could give the opposite 
outcome in the finger competition, that is, make the 'los- 
ing' finger for B = become the 'winning' finger when 
B > 0, for sufficiently generic sets of initial conditions. 

A satisfactory analytical understanding of the problem 
with B ^ has been achieved in two regimes: the ini- 
tial linear instability of the flat interface followed by the 
weakly non-linear regime [ p3| , and the asymptotic regime 
where surface tension selects the width of the single fin- 
ger |l5|, [l6| [It], The highly non- linear intermediate 
regime that connects the quasi-planar interface with the 
asymptotic single-finger regime has mostly been studied 
through numerical computation. The first exhaustive nu- 
merical studies were reported by Tryggvason and Aref 
p5| , who paid considerable attention to the influ- 
ence of viscosity contrast on the problem and studied 
both single-finger and multi-finger configurations. Later, 
Casademunt and Jasnow , 27 showed that the basin of 
attraction of the single-finger solution depends strongly 
on viscosity contrast and that only when one of the two 
fluid viscosities is negligible it can be claimed that the 
single finger is the universal attractor of the problem. In 
the present work we will restrict to this limiting case. 
DeGregoria and Schwartz p8l ^] observed that well- 
developed fingers split when surface tension is sufficiently 
decreased. This tip-splitting instability is related to the 
fact that the Saffman- Taylor finger is linearly stable but 
non-linearly unstable, and the size of the perturbation 
that triggers the tip-splitting decreases quickly with sur- 
face tension [|3(J. Dai and Shelley |3l| showed that for 
small B numerical computations are extremely sensitive 
to the precision used in the computations. As a conse- 
quence noise level has to be controlled with care in order 
to ensure that the computation is sufficiently accurate. 
Hou et al. j3^| developed a numerical method that deals 
with the numerical stiffness of the problem in an efficient 
manner, aiding the ability to perform long time computa- 
tions. More recently Ceniceros et al. 
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using very 

high precision arithmetic have been able to study the 
effect of extremely small surface tension in the circular 
geometry with suction, and they have observed that sur- 
face tension can produce complex ramified patters even 
without the presence of noise. 

An analytical treatment of this highly nonlinear and 
nonlocal free-boundary problem faces challenging diffi- 
culties. In particular, a perturbative study for small B is 
complicated by the ill-posedness of the zero surface ten- 
sion problem. Tanveer [ p"3"[ was able to overcome this ob- 
stacle by embedding the zero surface tension problem in 
a well-posed one. In addition, this well-posed extension 



of the B = problem allowed Baker et al. |35| to de- 
velop a numerical method to compute the time evolution 
of zero surface tension dynamics in a well-posed manner. 
Once the B = problem is formulated in a well-posed 
way the B ^ case can be studied using a perturbative 
approach. The main result of the asymptotic perturba- 
tive theory developed by Tanveer [|l3| is that the effect of 
surface tension may be manifest in a 0(1) time: the evo- 
lution of the same initial interface for B = and B ^ 
will in general differ after a time of order one, even if the 
B = solution is smooth for all time. Siegel et al. |2fJ| 
have extended the work of Ref. to later stages of the 
evolution, and through numerical computation for very 
small values of B they showed that smooth B = so- 
lutions are indeed significantly affected by the presence 
of arbitrarily small B in order-one time, thus confirming 
the predictions of the perturbative theory. The zero sur- 
face tension solutions studied by Siegel et al. |l^, 20 in 
the channel geometry were single-finger solutions with an 
asymptotic width A, specifically chosen to be incompat- 
ible with selection theory for vanishing surface tension. 
They found that the singular effect of surface tension was 
to widen the finger in order to reach the selected width. 
The surprising feature here is that the effect of surface 
tension is felt in order-one time, i.e., that the time lapse 
for which the regularized solution approaches the unper- 
turbed one as B — > is bounded. 

The present paper expands the work of Refs. [jlj^, |2(| 
in the spirit of Ref. |2^] , towards the study of multi-finger 
solutions. However, unlike the studies of jHl |2(J we chose 
zero surface tension multi-finger soltuions which are com- 
patible with asymptotic selection theory, that is, with an 
asymptotic finger width A = 1/2 — the selected value in 
the limit B — > 0. In this way we isolate the intrinsic fin- 
ger competition dynamics from the selection effects re- 
sponsible of restoring the asymptotic width. Two differ- 
ent kinds of two-finger zero surface tension solutions are 
studied, and in both cases it is shown that surface tension 
acts as a singular perturbation to the dynamics in order- 
one time, modifying dramatically the late time configu- 
ration of the interface not only quantitatively but also 
qualitatively. Specifically we show that paths in phase 
space associated with zero and nonzero surface tension 
evolution, and indeed the global topological structure of 
the the phase spaces, may differ appreciably, even for 
arbitrarily small B. In physical terms, our evidence sug- 
gests that the presence of arbitrarily small surface tension 
can completely alter outcome of finger competition when 
compared with the zero surface tension evolution. 

The paper is organized as follows: in Sec. II the equa- 
tions describing Hele-Shaw flow are introduced, and a 
class of two-finger zero surface tension solutions relevant 
to two-finger competition is presented and briefly dis- 
cussed. This class of solutions will be used as initial con- 
dition for numerical computation with B > 0. In Sec. Ill 
the basic features of the asymptotic theory are recalled, 
and the theory is applied to the zero surface tension solu- 
tions introduced in the previous section. The numerical 
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computations with finite (but small) B are presented in 
Sec. £v|. Sec. discusses and summarizes the results 
obtained in previous sections. 



II. ZERO SURFACE TENSION 

In this section we present the equations which gov- 
ern the interfacial dynamics in a rectilinear Hele-Shaw 
cell, following the formalism of (r| . We consider a class 
of exact, time-dependent zero surface tension solutions 
that are relevant to the finger competition problem, and 
briefly describe the solutions within this class. 

Consider Hele-Shaw flow in the channel geometry, in 
which a fluid of negligible viscosity displaces a viscous 
liquid. The equations governing the interfacial evolution 
can be conveniently formulated by first introducing a con- 
formal map z(C, t) which takes the interior of the unit 
semicircle in the C plane into the region occupied by the 
viscous fluid in the complex plane z — x + iy, in such a 
way that the arc C = e ls for s£ [0, w] is mapped to the 
interface and the diameter of the semi-circle is mapped to 
the channel walls (§3]. The mapping function z(£,t) has 
the form z((,t) = —(2/n)ln( + i + f((,t), and inside and 
on the unit semi-circle we require /(Cj*) to be analytic 
and zq(C,, t) ^ 0. In addition, we require that 



Im / = 



(1) 



on the real diameter of the semi-circle. This latter con- 
dition ensures that z maps the diameter to the channel 
walls. Under suitable assumptions (see the Schwartz 
reflection principle may be applied to show that / is an- 
alytic and Zf for \C\ < 1. 

The effective velocity field, averaged across the plate 
gap, is a two-dimensional potential flow satisfying 
Darcy's law 



V<p. 



(2) 



Here <p is a velocity potential defined by <p = — (6 2 /12 / u) p, 
where p is the pressure, /i is the viscosity and b is the 
gap width. Under the assumption of incompressibility 
(V • u = 0) the potential satisfies Laplace's equation 



V> = 0. 



(3) 



Incompressibility also implies the existence of a stream 
function tp. Therefore, one can define a complex velocity 
potential W{z, t) = <p + hp which is analytic for z in the 
fluid region of the channel. Its form as a function of C 
reads 



W(C,t) = -(2/n) lnC + i + co(C,i) 



(4) 



where w(C, t) is an analytic function inside the unit circle. 
The condition that there is no flow through the walls 
implies that Im u> — must hold on the real diameter 
of the unit semi-circle. In the absence of surface tension, 
u = (see Eq. (g)). 



At the interface we impose the usual boundary condi- 
tions. The kinematic boundary condition states that the 
normal component of fluid velocity at a point on the in- 
terface equals the normal velocity of the interface at that 
point, and takes the form 



Re 



Zl 



1 



■Re [CW C ] 



(5) 



The dynamic boundary condition specifies that the pres- 
sure jump across the interface is balanced by surface ten- 
sion, and is given by 



B 

Rew = — - — rRe 

kl 



1 + 



(0) 



The parameter B is the nondimensional surface tension 
and is defined by B = b 2 T/12fiVa 2 , where T is the sur- 
face tension, V is the fluid velocity at infinity and a is 
half the cell width. The equations given in (l)-(5) are in 
nondimensional form, with lengths and velocities nondi- 
mensionalized by a and V, respectively. 

When B = it is well known that pole singularities 
in zq (i.e. in fe) present in the exterior of the unit disk 
are preserved under the dynamics, i.e., such singularities 
are neither created nor destroyed, although the location 
of those which are initially present will evolve with time. 
Exact B = solutions consisting of a collection of pole 
singularities with constant amplitude have been the focus 
of extensive studies (see e.g. [|| ). The simplest such 
solution leading to nontrivial finger competition consists 
of a pair of singularites in the upper halfplane of |£| > 1, 
located at positions that are symmetric with respect to 
the y-axis. A second pair of poles conjugate to the first 
pair is required to satisfy the symmetry restriction (Q). 
This exact solution takes the form ^ 

z((, t) = -- InC + ~(1 - A + ie) ln(l - ) 

7T 7T (sit) 



) + d(t) + i (7) 



where A and e are real constants with < A < 1 and 
e > 0, and d(t) is real. The singularity locations are 
given by the complex parameter ( s (t), which satisfies a 
simple differential equation given in pl|. Analyticity of 
/(C)*) m the unit circle implies that |C*(*)| > 1- We 
employ the convention that Cs (t) is a complex number in 
the first quadrant. The amplitudes of the singularities, 
given here by the numbers 1 — A + ie and its conjugate, 
are chosen so that the asymptotic form of the solution 
consists of one or two steadily propagating fingers of total 
width A. The parameter e determines the nature of the 
finger competition for B = 0. 

The solution (Q) describes generically two different fin- 
gers of the nonviscous phase penetrating the viscous one. 
In the linear regime |Cs(*)| >> 1 the interface consists of 
a single bump or finger, and as time increases a second 
finger may develop and grow, depending on the value of 
ArgC s (0). 



4 



We summarize the features of the solution (|7|) that are 
most relevant to the study of finger competition. Con- 
sider first e = 0. In this case the asymptotic configuration 
consists of one or two fingers of total width A, depend- 
ing on the initial condition. The singularities move to- 
ward the unit disk, with the limit as t — > oo denoted by 
( s {t) — > e ld . When 9 — the asymptotic configuration 
is a single Saffman- Taylor finger growing in the center 
of the channel (this asymptotic configuration is denoted 
ST(R)), for 9 = tt/2 it is a 'side' Saffman-Taylor finger 
i.e. a pair of half fingers of total width A with tips lo- 
cated at the cell walls (denoted ST(L)), and for 9 = 7r/4 
it is a 'double' Saffman-Taylor finger, namely two identi- 
cal fingers of width A/2 with tips at x = 0, ±1 (denoted 
2ST). For any other value of 9 the asymptotic config- 
uration consists of two unequal steadily growing fingers. 
The two- finger asymptotic configuration is a consequence 
of the continuum of fixed points that is present in the 
phase portrait of the dynamical variables which govern 
the shape of the interface, namely (Re C(0 , Im£ s (t)). In 
order to correspond to the notation of pi] introduce the 
variable a(t) = a'(t)+ia"(t) = l/(i («(*))• Then the pla- 
nar interface corresponds to a = 0, the center Saffman- 
Taylor finger to a = —i, the side Saffman-Taylor finger 
to a = i, and the double Saffman-Taylor finger to a — 1 . 
Figure [l] shows the phase portrait of the dynamical sys- 
tem obtained from the substitution of the mapping de- 
fined by Eq. (@) into the evolution Eqs. (§,§) for B = 0, 
using the dynamical variables (a', a") . The asymptotic 
states with two advancing fingers corresponds in the dy- 
namical system (a', a") to a continuum of fixed points 
given by \a\ = 1. Therefore, for e = the solution (0) 
does not exhibit finger competition. In addition, it is im- 
portant to note that the evolution of (fjj) with e = is 
free of finite time singularities, i.e., ^ in the domain 
|C| < 1 for all time. 

For e ^ the continuum of fixed points is removed (see 
Fig. [I]), as is the double Saffman-Taylor finger fixed point. 
Consequently, the solution to Eq. (0) exhibits 'successful' 
competition, in the sense that the asymptotic interface 
shape consists of a single Saffman-Taylor finger or side 
Saffman-Taylor finger. The price to pay is the appear- 
ance of finite time singularities for a certain subset of 
initial conditions, in the form of a zero of impact- 
ing the unit disk (this is a generic feature of conformal 
map solutions zq composed of a finite number of pole 
singularities — see ||). Then, only the subset of initial 
conditions free of finite time singularities is capable of 
sustaining finger competition all the way to the t — > oo 
outcome. Nevertheless, one may ask whether the class 
of B = solutions that are free of finite time singu- 
larities may describe, at least qualitatively, the physical 
finger competition for positive surface tension in the limit 
B 0. 

In the following sections we focus on the class of initial 
data for which the B = solutions are devoid of finite 
time singularities, and examine the B > dynamics. We 
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FIG. I: (a) Phase portrait for A = 1/2 and e = 0. The region 
between the dotted lines corresponds to two-finger configura- 
tions and the other two regions to single-finger configurations, 
(b) Phase portrait for A = 1/2 and e = 0.1. For the points on 
the dashed line the two fingers have the same length. 

develop a general theory for how the presence of positive 
surface tension affects the outcome of finger competition. 
This will enable us to predict the winner of finger compe- 
tition, i.e., the eventual asymptotic state. Most interest- 
ingly, we find instances in which the presence of arbitrar- 
ily small surface tension leads to dramatically different 
outcomes in finger competition when compared with the 
zero surface tension evolution. 



III. ASYMPTOTIC THEORY 

Little is known about the effect of finite (but small) 
surface tension B on the dynamics of zero surface tension 
multifinger solutions, and in particular on the class of ex- 
act solutions (0). For single finger configurations, steady 
state selection theory predicts that the finger cannot have 
an arbitrary width. Indeed, for vanishing surface tension 
B — > the width A = 1/2 is selected, asymptotically in 
time. Thus, it is clear that surface tension has a criti- 
cal influence on single finger solutions with A =^ 1/2. The 
nature of this influence in the limit B — > has been inves- 
tigated by Siegel, Tanveer and Dai [|l9| ^0|, who present 
evidence that zero surface tension single finger solutions 
with A < 1/2 are significantly perturbed by the inclusion 
of an arbitrarily small amount of surface tension in order 
one time. The effect of surface tension is an increase of 
the finger width to reach the width predicted by selection 
theory. 

Consider now the effect of small surface tension on 
the exact (B = 0) two finger solution (^). When < 
B <C 1 the asymptotic perturbation theory developed in 
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Refs. |l9[ can be applied. This perturbation the- 
ory describes the effects of the introduction of a small 
amount of surface tension on initial data z((, 0) specified 
in the extended complex plane, i.e., in a domain including 
the 'unphysicaP region |£| > 1 (the extended domain is 
required to make the B — problem well-posed) . The ef- 
fect of finite B is most important near isolated zeros and 
singularities of z^(Q, 0), where a regular perturbation ex- 
pansion in B breaks down. (Away from these points the 
perturbation expansion is regular, at least initially.) For 
the class of solutions (Q) we are discussing, the isolated 
singularities of z^(C,, 0) are simple poles. The theory sug- 
gests that the introduction of finite surface tension mod- 
ifies the poles (£ s ) by transforming them into localized 
clusters of —4/3 singularities, but these clusters move at 
leading order according to the B = dynamics. Thus the 
effect of one of these clusters on the interface is approx- 
imately equivalent to that of the unperturbed (B = 0) 
pole-like singularity that has given birth to it. 

The influence of surface tension on the zeros of Zf(C, 0) 
is more complex. Each initial zero instantly gives birth to 
two localized inner regions, i.e., regions where the B = 
and B > solutions differ by 0(1). One of the two inner 
regions moves, at least initially, according to the B = 
dynamics of the original zero £o ]39| ] . Since the partic- 
ular zero surface tension solutions considered here have 
zeros that are either bounded away from the unit disk 
for all time or impact the unit disk only after long times, 
the inner region around Co(t) has a negligible influence 
on the interface. The second inner region created around 
Co(0) moves differently. The theory suggests that this in- 
ner region consists of a cluster of singularities, whose size 
scales like B 1 ^ 3 . Unlike the case discussed above this sec- 
ond inner region moves away from the B = zero since, 
to leading order in B, it moves like a singularity of the 
zero surface tension problem and this speed is different 
form the speed of the zero Co (t) which spawned the clus- 
ter. As this singularity cluster approaches the physical 
domain it may perturb the flow and the interface shape 



may differ significantly from that at B — shape. The lo- 
cation of this singularity cluster will be denoted by Cd(t), 
and following [[t3| we shall call it the daughter singular- 
ity. We emphasize that the dynamics of the daughter 
singularity cluster is determined at lowest order solely 
by the B = solution z°((, t), at least until it arrives at 
the surroundings of the unit circle, and therefore can be 
simply computed once the initial locations of the zeros of 
Z((£, 0) are determined. 

The daughter singularity evolution equation is given 
by (see @) 

&(«) = -<#(&(*),*); Ci(o) = Co(o) (8) 

where q® is defined by 



dCC- 


_ Re 


CW°(C',i)" 


c c 


-c 


* C °(C',*)I 2 



(9) 



and the superscript denotes that the function evalua- 
tions are done using the corresponding B = solution. 
The function — qi((, t) also gives the characteristic veloc- 
ity of a pole or branch point singularity of ^(C, t) located 
at position C in the region |C| > 1- The initial position 
£d(0) is a consequence of the fact that each zero Co(0) 
of the zero surface tension solution will give birth to a 
daughter singularity. From Eq. (||) it can be shown |ll| 
that d\(d\/dt < 0, so that the daughter singularity ap- 
proaches the unit circle and it can impact it in a finite 
time td, the daughter singularity impact time, satisfyng 
|Cd(irf)| = 1- In the limit B — > 0, the daughter singular- 
ity impact time td signals the time when the effects of 
the surface tension are felt on the physical interface. For 
times larger than td the B = interface and the B — > 
are expected to differ significantly. 

For the family of exact B = solutions the mapping 
function ([?]) has four pole-like singularities: ±C, and ±£ s , 
and four zeros ±Co+ and ±(o- of zq located at 



—2 



Co 2 + 



Co 2 - = 



-(A + *e)C 2 -(A- ie)C 
2(1 - 2A) 



-(A + ie )C s 2 -(A-ie)C s 



(\ + ie)a + (\-ie)C s +4(1-2A)|C.| 



2(1 -2A) 



-2 



(A + ie)Cf + (A-ie)< a +4(1 - 2A)|C| 



2(1 - 2A) 



2(1 - 2A) 



(10a) 
(10b) 



For the particular case A = 1/2 this solution presents 
only one pair of zeros ±£o located at 



C ( 2 



IC. 



2[(A + ze)C 2 + (A- l£ )C 



-r2 n 



(11) 



In the following it will be useful to define the real quantity 
f3 = — (A + ie)Cs — (A — Je)C s which appears in ( [To|) and 
(0)- 

Depending on the value of A the initial data may have 
zeros on both the real and imaginary axes, or all the 
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zeros may lie on a single axis. This difference has sig- 
nificant consequences in the finite surface tension dy- 
namics. M ore specifi cally , when A < 1/2 the zeros de- 
scribed in (10a) and ( 10b ) are located on both the real 
and imaginary axes of |£| > 1, namely at ±|£o+| and 
±z|£o-|- The situation is different for A > 1/2, which 
is further divided into two cases, depending on whether 
(3 2 + 4(1 - 2A)|C S | 2 > or < 0. In the former case all 
four singularities lie on the real axis (for (3 > 0) or on the 
imaginary axis (for (3 < 0). In the latter case the four 
zeros are located off the axes in conjugate pairs, i.e. at 
±£o and ±Co- Finally, when A = 1/2 the solution (0) has 
only two zeros, located on the real axis at ±\( s \ 2 /y/— 2/3 
when (3 < and on the imaginary axis at ±\£ s \ 2 /y/2/3 
when (3 > 0. Note that for A = 1/2 the B = solution 
has two less zeros than for A ^ 1/2. 

The initial zero locations described above have a criti- 
cal bearing on whether the daughter singularity will im- 
pact the unit disk flofl . Although all daughter singulari- 
ties approach the unit disk, their impact may be shielded 
by the presence of an inner region corresponding to a pole 
singularity. More precisely, since C,d and ( s obey the same 
dynamical equation, they will move together if they get 
close enough to each other. However, the inner region 
around a pole moves to leading order like the B = 
pole, i.e., it moves exponentially slowly toward |£| = 1 
when |£g| — 1 « 1, and does not impinge upon the unit 
disk in finite time 0. In this case the 0(B 1 ^ 3 ) inner re- 
gion around the daughter singularity will not affect the 
dynamics on |£| = 1, at least until t = 0(— Yn.B). Before 
this time, we expect the interface to be uninfluenced by 
the presence of the daughter singularity. This shielding 
mechanism is discussed in the context of single fingers in 
& 

Knowledge of the t — > oo asymptotic state and the ini- 
tial locations of zeros can be used to ascertain whether 
shielding can occur. The B = asymptotic state cor- 
responds to ( 2 (t — > oo) — > ±1. Thus, for A < 1/2, 
only one pair of daughter singularities may be shielded — 
never both — so at least one pair of daughter singular- 
ities will impinge on the unit disk. The daughter sin- 
gularities will also not be shielded when A > 1/2 and 
f3 2 + 4(1 - 2A)|C S | 2 < 0. However, for A > 1/2 and 
1 + 4(1 - 2A) |<"s 1 2 > it is possible for all the daugh- 
ter singularities to be shielded, since they lie on a single 
axis. The daughter singularities can also be completely 
shielded when A = 1/2. The different possibilities are 
schematically depicted in Fig. |^. 

We have numerically computed the daughter singular- 
ity impact time td for various values of A and e, using 
initial conditions close to the planar interface, \( s \ 2 = 20 
and various values of Arg[£ 2 ]. Figure ^ shows the phase 
portrait for different values of A and e with the daughter 
singularity impact indicated. From the plots it is imme- 
diately seen that for A < 1/2 at least one daughter sin- 
gularity always hits the unit circle, and for A > 1/2 some 
trajectories are free from daughter singularity impact. In 
addition, it is observed that for fixed A a larger value of 



(a) 



i C, (0) 




(b) 



-Cs(0) 



-Cs(0) 



Cs(0) 



Cs(0) 




FIG. 2: (a) Schematic representation of the dynamics of 
pole (£ s ) and daughter ((d) singularities for A < 1/2. (b) 
Schematic representation of one of the two possible dynamics 
of pole (£ s ) and daughter ((d) singularities for A > 1/2. 



e causes the daughter singularities to hit in shorter times 
(or less developed fingers) than a smaller value of e, and 
for fixed e larger A implies larger impact times. We have 
also checked that the daughter singularity impact occurs 
well before a finite time singularity, i.e., the impact of a 
zero of zq. Thus, the effect of surface tension is signifi- 
cant well before the curvature in the zero surface tension 
solution becomes large. 
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FIG. 3: Phase portraits for (a) A = 1/3 and e = 0.1, (b) 
A = 2/3 and e = 0.1, and (c) A = 1/3 and e = 1/2. The 
daughter singularity impact is indicated by the symbols. The 
+ symbol corresponds to the impact of (d+, X to the impact 
of (d- and * to the simultaneous impact of £d+ and fd_. 



It is noted that the A dependence of the daughter sin- 
gularity impact is consistent with the results of steady 
state selection theory [[151 16, [l?], [l8|]. According to se- 
lection theory, for small B the possible values of A are 
discretized: A must satisfy the relation A = A„(B), given 
to leading order by 



A„(B) 



1 + {\^C n B) 



2/3 



n = 0,l,2,. 



(12) 



where n parameterizes the branch of solutions. Note that 
A n > 1/2 for all n. The steady finger shape is to leading 
order a Saffman- Taylor finger, with the above values of 
A„ substituted for the width A. On the other hand for 
e > the asymptotic state of (Q) is a Saffman- Taylor 
finger of width A. From Eq. ( |l2| ) it is clear there exists 
a steady solution with width \ n {B) close to a Saffman- 
Tay lor finger of arbitrary width A > 1/2. Thus the shield- 
ing of the daughter singularity, which leads to the per- 
sistence of a Saffman- Taylor solution with A > 1/2 over 
long times, is consistent with steady state selection the- 
ory |ilj| . In contrast for A < 1/2 there are no nearby 
steady solutions. Thus, a Saffman Taylor finger with 
A < 1/2 cannot persist over a long time. We see that the 
impact of a daughter singularity provides a mechanism 
for the onset of finger competition, finger widening, and 
selection of a width A > 1/2. 

For e = the scenario is similar, except there is an 
added class of exact B > solutions. Magdaleno and 



Casademunt |36 have shown that two-finger solutions 
composed of steadily propagating but unequal fingers do 
exist for small nonzero B. The introduction of a small 
nonzero surface tension selects a discrete set of solutions 
from the continuum of fixed points of the B = phase 
portrait. The solutions are parameterized by the total 
width of the fingers A = Ai + A2 and the relative width 
q = Ai/A, and the introduction of finite B discretizes 
the possible values of the parameters. In particular, they 



must satisfy a condition of the form A = \ n {B) and q = 
Qn,m(B) where n and m are integers. The expression for 
\ n (B) at lowest order is equivalent to Eq. (|l2|), but with 
different coefficients C n . The shape of these solutions 
are given to leading order (in the limit t — > 00) by (j?]) 
with allowed value of X n {B) substituted for the width 
A. Again, \ n {B) > 1/2, and the consistency between 
daughter singularity impacts and steady state selection 
theory follows as above. 

We conjecture that the outcome of interfacial shape 
evolution after the daughter singularity impinges is in 
general independent of the particular finger on which 
the impact first occurs i.e., independent of the point at 
which Q(t) impacts on |£| = 1. More specifically, we sur- 
mise that impact on either the shorter (trailing) or larger 
(leading) finger retards the velocity of that finger, and is 
accompanied by the widening of the leading finger, so as 
to maintain a constant fluid flux at infinity. The widened 
leading finger then shields the trailing finger, preventing 
it from further growth. Thus, the finger which is leading 
at the time of the daughter singularity impact 'wins' the 
competition, in the sense that it will evolve for t — > 00 
to the ST steady finger. To examine this conjecture and 
study the dynamics of finger competition with finite (but 
small) surface tension we have numerically computed the 
evolution of an interface with initial conditions given by 
the conformal mapping Eq. (Q) close to the planar inter- 
face (|Cs(0)|" 2 <C 1). The results are reported in the next 
section. 



IV. NUMERICAL RESULTS 

Numerical computations have been performed for B > 
0, using an initial interface corresponding to the explicit 
B = solutions discussed in Sec. II. The effect of positive 
surface tension on this class of solutions is explored for 
various values of e and a variety of initial pole positions. 

We employ the numerical method introduced by Hou et 
al. |32| and used in other studies of small surface tension 
effects in Hele-Shaw flow g|j20|, || |§. The method is 
described in detail in Ref. |32|. It is a boundary integral 
method in which the interface is parameterized at equally 
spaced points by means of an equal-arclength variable a. 
Thus, if s{a,t) measures arclength along the interface 
then the quantity s a (a,t) is independent of a and de- 
pends only on time. The interface is described using the 
tangent angle 9(a,t) and the interface length L(t), and 
these are the dynamical variables instead of the interface 
x and y positions. The evolution equations are written in 
terms of 9(a, t) and L(t) in such a way that the high-order 
terms, which are responsible of the numerical stiffness of 
the equations, appear linearly and with constant coef- 
ficients. This fact is exploited in the construction of an 
efficient numerical method, i.e., one that has no time step 
constraint associated with the surface tension term yet is 
explicit in Fourier space. We have used a linear propaga- 
tor method that is second order in time, combined with a 
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spectrally accurate spatial discretization. Results in this 
section are specified in terms of the scaled variables 

t = irt; B — tt 2 B\ x — nx; y — ny. (13) 

instead of the original ones used in previous sections. 

The number of discretization points is chosen so that 
all Fourier modes of 8(a,t) with amplitude greater than 
round-off are well resolved, and as soon as the ampli- 
tude of the highest-wavenumber mode becomes larger 
than the filter level the number of modes is increased, 
with the amplitude of the additional modes initially set 
to zero. The time step At is decreased until an additional 
decrease does not change the solution to plotting accu- 
racy, nor lead to any significant differences in any quan- 
tities of interest. In a typical calculation 512 discretiza- 
tion points are initially used, and the initial time step is 
At = 5 • 10~ 4 . For small values of surface tension numer- 
ical noise is a major problem, and the spurious growth of 
short-wavelength modes induced by round-off error must 
be controlled. To help prevent this noise-induced growth 
at short wavelengths spectral filtering J3?j is applied. Ad- 
ditionally, we minimize noise effects and also assess the 
time at which these effects become prevalent by employ- 
ing extended precision calculations, as described in the 
next section. 

Our main interest is to uncover the role of surface ten- 
sion in the dynamics of finger competition. To isolate 
the features of finger competition from those of width 
selection, we will concentrate on B = solutions with 
A = 1/2, the value selected by surface tension in the 
limit B — > 0. Since the B = dynamics for e = and 
e 7^ is quite different the numerical results for the two 
cases will be presented separately. 

A. Solutions with e — 

We first consider parameter values X — 1/2 and e = 0. 
A typical set of interfacial profiles is shown in Fig.0. The 
initial data is given by the mapping function Eq. (H), with 
A = 1/2, e = 0, d(0) = and ( s 2 (0) = 20 exp(i7r/6). With 
this value of Cs(0) the initial interface is well inside the 
linear regime. Evolutions are shown for different values 
of B, and the B = interface evolution is also plotted 
for comparison. In all these evolutions the filter level is 
set to 10~ 13 , although later we shall make comparisons 
to profiles computed at higher precision. 

For the largest value of surface tension the computed 
B > and the exact B — solutions first differ appre- 
ciably at the seventh curve, corresponding to t 3. At 
this point the velocity of the small finger (at the chan- 
nel sides) begins to decrease and it is clearly left behind 
when compared with the small finger evolution in the 
B = solution. Eventually, the advance of the small fin- 
ger is completely suppressed and the larger finger widens 
to attain a width close to 1/2 of the channel. For a 
smaller value of surface tension, for instance B = 0.001, 
the evolution displays qualitatively the same behavior. 



The B > interface differs appreciably from the B = 
sightly later than before (i.e., at the eighth curve) and 
the region where the two solutions differ most is to some 
extent more localized around the small finger than for 
larger values of B. Additionally, for this value of surface 
tension the effect of numerical noise is clearly exhibited 
in the interfacial profiles. Here the tip-splitting and side- 
branching activities are a clear effect of numerical noise, 
as can be easily checked redoing the computation with a 
different noise filter level. 

In order to suppress or delay the branching induced by 
numerical noise that appears for small values of surface 
tension it is necessary to use higher precision arithmetic, 
e.g. quadruple precision (128- bit arithmetic). The fil- 
ter level can then be reduced by a large amount and the 
outcome of spurious oscillations is substantially delayed. 
Figure || shows the effect of reducing the filter level to 
10~ 27 . The B = solution is plotted, as well as the 
computation with double precision. For B = 0.001 the 
branching is totally suppressed, at least for the times we 
have computed, but for smaller values of B the use of 
quadruple precision is only able to delay the branching 
and not totally suppress it. The quadruple precision com- 
putation confirms the results observed with lower preci- 
sion: the introduction of finite (but small) surface ten- 
sion results in the suppression of the small finger. From 
Fig. H one can also see that for long times, when the 
interface is clearly affected by numerical noise (in the 
double precision curve), the noise- induced branching is 
restricted to the large finger, and the small finger is basi- 
cally unaffected by noise. This observation suggests that 
the small finger shape, as well as its tip velocity and tip 
curvature, can be trusted even when the large finger has 
developed tip-splittings and side-branchings due to the 
spurious growth of round-off error. 

Figure ^ shows the tip velocity of both fingers ver- 
sus t for decreasing values of surface tension. It can be 
seen that the velocity of the large finger is only slightly 
affected by surface tension, whereas the velocity of the 
small finger is substantially reduced by the inclusion of 
finite B. As B is decreased the tip velocity of the small 
finger is more faithful to the B = evolution before 
the daughter singularity impact (shown by a cross), and 
clearly veers away from the B = velocity later in the 
evolution, consistent with asymptotic theory. Note that 
at the smallest value of B the tip velocity of the large 
finger drastically differs from the B = velocity at late 
times. This discrepancy is a manifestation of noise effects 
in the neighborhood of the large finger tip. However, as 
previously seen, the small finger is basically unaffected 
by noise at the times we have plotted. 

In order to further verify that the daughter singular- 
ity impact is responsible for the observed change in the 
small finger tip speed we follow the scheme introduced in 
|p0f . Define t p as the time when the computed tip veloc- 
ity differs by p from the B = tip velocity. According to 
asymptotic theory this t p will be a linear function of B 1 / 3 
in the limit B — > as long as p is small enough. Figure R 
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(a) (c) 




FIG. 4: Evolution of an initial condition of the form (Q) with A = 1/2, e = and Cs(0) = 20 exp(i7r/6). The solid lines 
correspond to surface tension B values (a) 0.01, (b) 0.005, (c) 0.001 and (d) 0.0005. The dashed lines correspond to the zero 
surface tension evolution. The time difference between different curves is 0.5. The physical channel in the y direction extends 
from the origin to the dotted line, and the region above is plotted for better visualization of the lateral finger. 



shows t p versus £}V3 for various values of p, and it can 
be seen that t p exhibits the predicted behavior. More- 
over, we have extrapolated the B = value of t p using 
the two points of lowest B and the result is very close to 
id, whose value is represented by a cross. We conclude 
that the impact of the daughter singularity is associated 
with the dramatic change of the B > solution when 
compared to the zero surface tension solution, reducing 
the velocity of the small finger and eventually suppress- 
ing it. In contrast, for the B = dynamics the small 
finger 'survives', propagating with the same asymptotic 
speed as the larger finger. Note that the average inter- 
face advances at unit velocity, and a tip velocity below 1 
implies that the finger is retreating in the reference frame 
of the average interface. 

In summary then, our numerical results show that the 
computed interface for B ^ follows the B — evolution 
for an O(l) time interval — roughly corresponding to the 
daughter singularity impact time — and that at further 
times the velocity of the small finger decreases while the 
large finger widens. The small finger eventually comes to 
a halt and the larger (leading) finger reaches an asymp- 
totic width slightly above 1/2, the width singled out by 
selection theory. It is noted that for the initial condition 
we have studied the daughter singularity impact takes 
place on the tip of the small finger. Therefore, the in- 
fluence of surface tension on the interface should be sig- 



nificant first around the impact point, that is, the small 
finger tip. Our numerical results show that in fact this 
is the case; the initial effect of the daughter singularity 
impact is to slow and then completely stop the growth 
of the small finger. Later on, as the singularity cluster 
centered in Qd spreads over the unit circle, the effect of 
surface tension is felt by the whole interface and the large 
finger widens to reach the selected width. 

We have also studied the finite surface tension dynam- 
ics for a more general class of initial conditions. More 
precisely, we have studied initial conditions of the form 
C S 2 (U) = 20exp(i rm/12) where n = 0,±1,...,±6, and 
have obtained the same qualitative results as in the case 
previously studied, namely that the presence of small sur- 
face tension suppresses the growth of the finger which is 
trailing at the time of daughter singularity impact. In 
order to compare the B = and the B ^ dynamics 
in a compact and global way we have plotted the phase 
portrait for B — using the the tip velocities v\, V2 as 
dynamical variables. In the laboratory frame they read 

i + KC 2 -C 2 ) + C 2 C 2 hA , 

Vl = a 2 + i(c 2 -c 2 )/2 (14a) 

i - KC 2 - C 2 ) + c 2 c 2 nAU , 

V2= c 2 c 2 -i(c 2 -c 2 )/2 - (14b) 

Now a comparison between dynamics for B = and 
B 7^ is straightforward since the trajectories can be 
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FIG. 5: Evolution of an initial condition of the form ([]) 
with A = 1/2, e = and < s 2 (0) = 20 exp(i7r/6). The solid line 
correspond to B — 0.001 with a filter level equal to 10 -27 , the 
dotted line corresponds to the same B but with the filter level 
equal to 10~ 13 and the dashed line corresponds to the zero 
surface tension solution. The time difference between curves 
is 0.5. As in Fig. ^, the physical channel in the y direction 
extends from the origin to the dotted line. 
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FIG. 6: Computed tip velocities for the initial condition of 
Fig. |H| , (a) corresponds to the central (large) finger and (b) 
to the lateral (small) finger. The daughter singularity impact 
time td is indicated by the + symbol. The value of B is: 
(solid line), 0.0002 (dotted line), 0.0005 (dashed line), 0.001 
(long dashed line), 0.005 (dot-dashed line) and 0.01 (dot-dot- 
dahed line). 
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FIG. 7: The time t p (defined in the text) versus B 1/3 . From 
top to bottom, p=0.1, 0.05, 0.02, 0.01, 0.005. The daughter 
singularity impact time td is inducated by a x symbol, and 
the curves are linearly extrapolated for comparison. 



plotted together and compared. In addition, the tip ve- 
locity is a useful variable because it contains geometric 
information; specifically the inverse of the tip velocity 
is equal to the width of the finger in the asymptotic 
(t — > oo) regime. It is important to note that (1)1,1)2) 
are dynamical variables for the B = problem, so that 
the plot of the zero surface tension trajectories onto the 
space (v±, V2) is a true phase portrait. On the other hand 
(ui, V2) are not state variables of the problem with finite 
surface tension, so in this case we simply obtain a projec- 
tion onto the (v\,V2) space of the original B 7^ trajec- 
tory, which is embeded in the infinite-dimensional phase 
space of interface configurations. 

Figure || shows the phase portrait for B = together 
with the tip velocities obtained from the initial condi- 
tions described above for B = 0.01. From the figure it 
is evident that the introduction of finite surface tension 
has substantially changed the global phase dynamics of 
the problem. Only one B = 0.01 trajectory connects the 
planar interface (1,1) and the 2ST point (2,2), corre- 
sponding to the unsteady double Saffman- Taylor finger. 
Any other B = 0.01 trajectory ends in one of the two 
ST finger points, ST(L) at (2,0) and ST(R) at (0,2). In 
contrast, the (2, 2) point, equivalent to the continuum 
of fixed points present with the (a 1 , a") or (Re^ s ,Im^ s ) 
variables, has a finite basin of attraction for B = 0. The 
introduction of finite surface tension has dramatically 
changed the zero surface tension (1)1,1)2) trajectories, to 
the extent that the B = phase portrait and the B ^ 
projection are not topologically equivalent. This result is 
not a complete surprise, since it was anticipated from the 
structural instability of the dynamical system governing 
the evolution of Eq. (0) for e = ]2lJ . A more dramatic 
example of topological inequivalence of phase portraits 
will be given in the next subsection, when we consider 
the case e 7^ 0. 
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FIG. 8: Plot of the evolution of initial conditions of the 
form (@) with A = 1/2, e = and C s 2 (0) = 20exp(i nvr/12) 
and n = 0, ±1, ±6 in the (vi, V2) or tip speed space. The 
solid line corresponds to B = 0.01 and the dashed line to 
B = 0. 




a' 



FIG. 9: Comparison between the B = trajectories and 
the projected evolutions with B = 0.01, for the initial con- 
ditions of Fig. [|. The solid line corresponds to B = and 
the dashed lines to the projection of the B = 0.01 evolutions. 
The daughter singularity impacts are indicated by a circle. 



Although the use of the variables (ui,i>2) has allowed 
us to project the finite surface dynamics onto the zero 
surface tension phase portrait this projection has one 
major limitation: it only considers a local quantity, the 
tip velocity. We have also considered a projection that 
takes more global properties of the interface into account. 
Specifically, given a computed B ^ solution for an ini- 
tial condition of the form (0), one can use a suitable 
norm to define a 'distance' between the computed inter- 
face and the B — interface obtained from the mapping 
function Eq. ([?]). We choose this 'distance' to be the 
area enclosed between the two interfaces at a given time. 
Additionally, we define a projection of the B ^ in- 
terface onto the B = phase space (with phase space 
variables (ReCs,ImCs)) by selecting the value of £ s that 
minimizes the 'distance' between the two interfaces, with 
the restriction that the position of the two mean inter- 
faces must be the same. The latter condition ensures 
that the projection satisfies mass conservation. 

Figure || shows the B — phase portrait and the cor- 
responding projected evolution for surface tension B = 
0.01. Again, the plot clearly shows that the introduc- 
tion of finite surface tension modifies the phase portrait 
of B = 0. The projected trajectories are initially close 
to the B = dynamics, but for well developed fingers 
(corresponding to |a| ~ 1) the projection departs from 
the B = trajectory towards the Saffman- Taylor fixed 
point, located at a' — Q,a" = 1. The projected tra- 
jectory only remains close to the corresponding B = 
trajectory when the latter evolves towards the Saffman- 
Taylor fixed point. More precisely, the continuum of fixed 
points present for B — has been removed by surface ten- 
sion and the Saffman- Taylor fixed point is the universal 
attractor of the dynamics for finite surface tension. 

In Fig. [l(] the projection for decreasing values of B is 
plotted, using the initial condition C|(0) = 20 exp(i7r/6). 
As B is decreased the projected trajectory gets closer 
to the B = trajectory, but as it approaches the point 
when the daughter singularity impinges the unit circle 
(this point is signaled by a cross) the projection departs 
from the B = trajectory and approaches the Saffman- 
Taylor fixed point, consistent with asymptotic theory. 



B. Solutions with e / 

The continuum of fixed points present for e = is 
absent for e 7^ 0, but in this case finite-time singular- 
ities in the form of zeros of zq impinging on the unit 
disk do appear for some initial conditions. Therefore, we 
can expect that the effect of finite surface tension will be 
somewhat different than for e = 0. Firstly, the presence 
of surface tension should eliminate finite-time singular- 
ities, and secondly, finite B could modify the basin of 
attraction for the two attractors of the B = dynamical 
system, namely the side Saffman- Taylor finger and the 
center Saffman- Taylor finger. 

To explore this, we have performed computations with 
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FIG. 10: Comparison between the B = trajectories and 
the projection of the evolution of the initial condition given 
by Eq. (§) with A = 1/2, e = and C*(0) = 20 exp(i7r/6), 
where ▲ corresponds to B = 0.001, o to B — 0.005, □ to 
B — 0.01 and x to B — 0. The daughter singularity impacts 
are indicated by a plus. 
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FIG. 11: Plot of the evolution of initial conditions of the 
form (|) with A = 1/2, e = 0.1 and Cf (0) = 20exp(i nvr/12) 
and n — 0, ±1, ±6 in the (vi,V2) or tip speed space. The 
solid line corresponds to B — 0.01 and the dashed line to 
B — 0. The computed trajectory that most nearly separates 
the two basins of atraction is also plotted. Note that the 
long time behavior of the third and fourth B = 0.01 curves 
(counting from the upper left trajectory in clockwise direc- 
tion) is dramatically different from the corresponding B = 
solutions. 



A = 1/2 and e = 0.1 with initial conditions Cs(0) = 
20exp(i nn/12) and n — 0,±1,...,±6. Initially we set 
B = 0.01 and use a value of the noise filter level equal to 
10~ 13 , which suffices due to the relatively large value of 
B. The easiest way to compare both dynamics, finite B 
and B = 0, i s to plot their trajectories in velocity space. 
Thus, in Fig. 
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the tip velocities (i>i,U2) of the B = 0.01 
computation are plotted together with the tip velocities 
for B = 0. For arbitrary e and A the tip velocities of the 
B = solution read 



L'2 



i + KC 2 -C 2 ) + C 2 C 



2X2 



C 2 C 2 



<C 2 + C 2 ) + iA(C 2 -C 2 ) 
i-i(C 2 -C 2 ) + C 2 C 2 



(1-2A) 



CICI + e(C 2 + C 2 ) - iA(C 2 - C 2 ) - (1 - 2A) 



(15a) 



(15b) 



From the plot one can see that most B — 0.01 velocity 
trajectories follow (at least qualitatively) their B = 
counterparts, in the sense that they end up in the same 
fixed point. However, the second, third and fourth trajec- 
tories (counting from the upper left trajectory in clock- 
wise direction) differ significantly from their B = coun- 
terparts. The second B = 0.01 trajectory moves apart 
from the B — solution simply because the latter de- 
velops a finite-time singularity, which is regularized by 
the introduction of finite surface tension. However, the 
third and fourth trajectories exhibit a quite surprising 
behavior: the computed interface with B = 0.01 ends 
up in a different fixed point than the exact B = solu- 
tion, despite the fact that the B = solution is smooth 
for all time and has the asymptotic width that would be 
selected by vanishing surface tension. 

In order to get further insight into this behavior 
we have computed the evolution for decreasing values 
of B using the specific initial pole position C 2 (0) = 
20exp(— i7r/6), with A = 1/2 and e = 0.1. Quadruple 
precision has been used when it has been necessary. Fig- 
ure [l2] shows its evolution for four values of the surface 
tension parameter, together with the B = solution. 
The differences between the two interfaces for long times 
are readily apparent. When B = the finger in the cen- 
tral position stops growing and the side finger wins the 
competition, whereas for B > we encounter the oppo- 
site situation — namely, the central finger surpasses the 
side finger and wins the competition. For the smaller 
values of B the finger on the sides has not quite stopped 
growing when the computation is stopped, although its 
tip speed shows a marked decrease over that for B = 
and is less than that of the central finger. The side fin- 
ger tip speed is also decreasing at the final stage of the 
computation. The tip speed trend in the limit B — > 
is further illustrated in Fig. (|l3|). This figure shows the 
tip speed versus time of each finger for a sequence of de- 
creasing B. The plot suggests that upon impact of the 
daughter singularity the side finger velocity levels off and 
eventually decreases, whereas the velocity of the center 
finger is nearly unaffected and continues to increase. The 
trend is indicative of the center finger "winning" the com- 
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FIG. 12: Evolution of an initial condition of the form ((?] ) with A = 1/2, e = 0.1 and C s 2 (0) = 20 exp(-i7r/6). The solid lines 
correspond to surface tension B values (a) 0.01, (b) 0.005, (c) 0.001 and (d) 0.0005. The dashed lines correspond to the zero 
surface tension evolution. The time difference between different curves is 0.5. The physical channel in the y direction extends 
from the origin to the dotted line. 



petition in the B > dynamics, while the opposite occurs 
for B = 0. Finally, it is noted that the influence of sur- 
face tension is first felt by the smaller finger, which is the 
recipient of the daughter singularity impact. Afterwards 
the leading finger begins to widen, in a manner consis- 
tent with the conjecture in Sec. III. Further remarks on 
this point are made in Sec. V. 

The projection method described in the previous sec- 
tion has been also applied to this case, and the results 
are displayed in Fig. [14] in the particular case B = 0.01. 
It can be seen that for most trajectories the projection 
stays close to the B — curves, even for long times. The 
daughter singularity impact still leads to O(l) differences 
between the B = and B > solutions, although the 
impact does not produce changes in the outcome of finger 
competition. However, as expected some of the trajecto- 
ries (namely the third and fourth as measured clockwise 
from the bottom) do indicate significant qualitative dif- 
ferences in the long time evolution. The plot provides a 
simple depiction of the topological inequivalence of the 
B > and B — dynamics [ ft2| . 

It has been shown that the introduction of a finite B 
has not changed the attractors of the problem, but it has 
changed their basins of attraction. Interestingly, in the 
B = case there does not exist a single separatrix tra- 
jectory between the two Saffman- Taylor attractors, but 
rather a finite region, corresponding to the set of trajec- 



tories ending in cusps, that acts as an effective separatrix. 
Since for finite surface tension there are no cusps, it can 
be assumed that there is a single trajectory that separates 
the two basins of attraction. Obviously, this trajectory 
will depend on the value of the surface tension parameter. 
More precisely, the initial condition Cs (0) corresponding 
to the separatrix trajectory will be a function of the sur- 
face tension B. To quantitatively characterize this set 
of initial conditions we have studied the dependence of 
the separatrix trajectory in a neighborhood of the planar 
interface fixed point as a function of B, using initial con- 
ditions of the form £f (0) = 20 exp(i9) . For a given initial 
condition £ s (0) introduce the parameter 9 sep (B), defined 
as the unique value for which the evolution is attracted 
toward the fixed point ST(L) when 9 > 9 sep and to the 
fixed point ST(R) when 9 < 9 sep . 

Figure [l5] shows the plot of 9 sep versus B, and it is 
observed that as B decreases, 9 sep saturates to a fixed 
value, namely 9 sep (B -> 0) = -0.4843 ± 0.0009. It is 
interesting to compare this value to the position of the 
separatrix region for B — 0, which is located between 
6f=° = -0.95758 and 9 B=0 = -1.04796. The separatrix 
for finite B lays outside and far away from the separatrix 
region for B = 0, even for vanishing surface tension. Our 
evidence therefore suggests that any B = trajectory 
located between the trajectories defined by 9 sep (B — > 0) 
and 9 B=0 will not describe, even qualitatively, the reg- 



14 




(b) 



2.5 




FIG. 13: Computed tip velocities for the initial condition 
of Fig. [l2] , (a) corresponds to the central finger and (b) to 
the lateral finger. The daughter singularity impact time td is 
indicated by the + symbol. The value of B is: (solid line), 
0.0002 (dotted line), 0.0005 (dashed line), 0.001 (long dashed 
line), 0.005 (dot-dashed line) and 0.01 (dot-dot-dashed line). 
The deviations observed at late times for B — 0.0002 and 
B — 0.0005 in (b) are due to numerical noise. 



ularized dynamics in the limit B — > 0, since the finger 
that will 'win' the competition under the B = dynam- 
ics will 'lose' under the B — ► dynamics. Thus, there 
exists a positive measure set of initial conditions of the 
form (^) such that the evolution with B — ► cannot be 
qualitatively described by its evolution under B = dy- 
namics. This is a dramatic consequence of the singular 
nature of surface tension on the dynamics of finger com- 
petition which is not related to steady state selection, but 
confirms the ideas of the proposed dynamical solvability 
scenario in Ref. [E2|. 



SUMMARY AND CONCLUDING REMARKS 




FIG. 14: Comparison between the B = trajectories and 
the projected evolutions with B = 0.01, for the initial condi- 
tions of Fig. [l2|. The solid line corresponds to B — and the 
dashed lines to the projection of the B — 0.01 evolutions. The 
daughter singularity impacts are indicated by a circle. Note 
that the fourth R ~> tra.iectorv fas measured counterclock- 
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The asymptotic theory developed in Refs. (f| |(J pre- form (0) with A = 1/2, e = 0.1 and C?(0) = 20exp(i 
diets the existence of regions of the complex plane where 
the zero surface tension solution and the finite surface 
tension solution differ by 0(1). These regions are the 
daughter singularity clusters, and their influence is felt in 
the physical interface when they are close to the unit cir- 
cle. Daughter singularities move towards the unit circle, 
and when their motion is not impeded by other singular- 
ities they reach the unit circle in O(l) time. When the 



distance between the daughter singularity and the unit 
circle is 0(B 1 ^ 3 ) the interface can display 0(1) discrep- 
ancies with respect the interface of the B — solutions. 
However, the asymptotic theory does not predict the na- 
ture of the discrepancies caused by daughter singularity 
impact. 
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Siegel et al. [^0| showed numerically that the effect of 
the daughter singularity impact on the tip (in a single- 
finger configuration with A < 1/2) was a retard in 
the velocity of the finger accompanied by a widening 
of the finger. However, this provided small insight on 
the effect of the impact in multi-finger configurations, 
where finger competition could be substantially affected 
by the presence of finite surface tension, as suggested in 
Refs. g^g. 

Since the precise effect of the daughter singularity can- 
not be established by the asymptotic theory it is neces- 
sary to use numerical computation in order to establish 
the effects of daughter singularity on the dynamics of 
the interface. We have focused our efforts on uncovering 
the role of surface tension in the dynamics of two finger 
configurations, which is the simplest situation exhibiting 
nontrivial finger competition. Two different types of two- 
finger zero surface tension solutions have been studied. 
The first type (e = 0) does not exhibit finger competition 
when B = but rather contains asymptotic configura- 
tions consisting of two unequal steady fingers advancing 
with the same speed. These two-finger steady state so- 
lutions form a continuum of fixed points in the phase 
space of the corresponding, reduced dynamical system, 
which is structurally unstable. Numerical computations 
with small surface tension show that the introduction of 
a small B removes the continuum of fixed points and 
triggers the competition process which was absent for 
B = by restoring the saddle-point (hyperbolic) struc- 
ture of the appropriate multifinger fixed point. The sec- 
ond type (e ^ 0) of two-finger solution we have studied 
exhibits finger competition for B = 0, but the numerical 
computation with small B has shown that the long time 
configuration of the computed interface may be quali- 
tatively different from the B = solution for a broad 
set of initial conditions, in the sense that the finger that 
'wins' the competition is not the same with and without 
surface tension. Thus, the presence of surface tension 
seemingly can change the outcome of finger competition 
even in configurations that are well behaved and smooth 
for all time and whose asymptotic width is fully compat- 
ible with the predictions of selection theory for vanishing 
surface tension. This unexpected result shows that sur- 
face tension is not only necessary to select the asymptotic 



width and to prevent cusp formation, but plays also an 
essential role in multifinger dynamics through a drastic 
reconfiguration of the phase space flow structure. 

Our calculations support the conjecture that impact 
on either the shorter or larger finger retards the velocity 
of that finger, and is accompanied by the widening of the 
larger finger. As a consequence, in general the outcome 
of finger competition is independent of the particular fin- 
ger on which the impact first occurs, and the finger which 
is leading at the time of the daughter singularity impact 
'wins' the competition. This recipe fails only for inter- 
facial configurations with very similar fingers, when not 
only the position of the finger (which finger is leading) 
but also the tip velocities (a trailing finger can have for 
a certain time a larger velocity than the leading one) at 
the impact time may play a role. 

The main conclusion of the present work is that surface 
tension is essential to describe multifinger dynamics and 
finger competition, even when the corresponding zero 
surface tension evolution is well behaved and compatible 
with selection theory. That is, we have detected singular 
effects of surface tension on the dynamics of finger 
competition that are not directly related to steady state 
selection. These can be properly interpreted in the 
context of an extended dynamical selection scenario as 
described in Ref. |2^] where the reconfiguration of phase 
space flow by surface tension can be traced back to the 
restoring of hyperbolicity of multifinger fixed points. 
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